import numpy as np
from numpy import linalg

class PCA:

    '''
    dataset 形如array([样本1,样本2,...,样本m]),每个样本是一个n维的ndarray
    '''
    def __init__(self, dataset):
        self.indata = np.matrix(dataset, dtype='float64')
        self.ceprei_const_axis = 0 # 计算每列或每行的均值，0：计算每列均值，1：计算每行均值
        self.ceprei_const_rowvar = False # 是否每一行代表一个特征，每一列代表一个样本

    '''
    求主成分;
    threshold可选参数表示方差累计达到threshold后就不再取后面的特征向量.
    '''
    def principal_comps(self, threshold = 0.85):

		# 标准化
        X_menu = self.indata - np.mean(self.indata, axis=self.ceprei_const_axis)
        # 求 协方差矩阵
        Cov = np.cov(X_menu, rowvar=self.ceprei_const_rowvar)
		# 求特征值和特征向量
        eigValues, vectors = linalg.eig(Cov)
        # 按照特征值从大到小排序
        sorted_indices = np.argsort(eigValues)[::-1]
        eigValues = eigValues[sorted_indices]
        vectors = np.round(vectors[:, sorted_indices], 4)
        # 计算方差贡献率
        var_exp = np.round(eigValues / np.sum(eigValues), 4)
        # 计算选取的主成分
        proj_matrix = vectors
        # 降维, 根据累积贡献率来选取
        cum_var_exp = 0
        for i, var in enumerate(var_exp):
            cum_var_exp += var
            if cum_var_exp >= threshold:
                proj_matrix = proj_matrix[:, :i + 1]
                break

        X_pca = np.dot(self.indata, proj_matrix)
        # 纠正主成分正负号问题
        for i in range(proj_matrix.shape[1]):
            if np.sum(proj_matrix[:, i]) > 0:
                proj_matrix[:, i] = np.round(-proj_matrix[:, i], 4)
                X_pca[:, i] = np.round(-X_pca[:, i], 4)

        return proj_matrix, var_exp, X_pca, vectors

if __name__ =="__main__":
    # test_pca_with_txt_data();
    dataset = [
         [66, 64, 65, 65, 65],
         [65, 63, 63, 65, 64],
         [57, 58, 63, 59, 66],
         [67, 69, 65, 68, 64],
         [61, 61, 62, 62, 63],
         [64, 65, 63, 63, 63],
         [64, 63, 63, 63, 64],
         [63, 63, 63, 63, 63],
         [65, 64, 65, 66, 64],
         [67, 69, 69, 68, 67],
         [62, 63, 65, 64, 64],
         [68, 67, 65, 67, 65],
         [65, 65, 66, 65, 64],
         [62, 63, 64, 62, 66],
         [64, 66, 66, 65, 67]]
    p = PCA(dataset)

    proj_matrix, var_exp, X_pca, vectors = p.principal_comps(threshold=0.95)
    print("特征向量:\n", vectors) # 特征向量
    print("计算选取的主成分txt:\n", proj_matrix) # 计算选取的主成分
    print("计算方差贡献率:\n", var_exp) # 计算方差贡献率
    print("降维后的数据:\n", X_pca) # 降维
  #  print(lst)
